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PSEUDOSPECTRAL  SOLUTION  OF  INVISCID  FLOWS 
WITH  MULTIPLE  DISCONTINUITIES 


I.  INTRODUCTION 

/ 

The  author  has  shown  [1-3]  that  a  pseudospectral  technique  may  be 
coupled  with  fourth-order  artifical  viscosity  and  spectral  filtering  [4]  to 
solve  inviscid  flow  fields  in  which  a  single  discontinuity  is  present.  The 
flow  fields  treated  in  this  manner  have  been  both  one  and  two  dimensional  in 
character;  the  former  consisting  of  a  shock  wave  propagating  in  the  co¬ 
ordinate  direction  and  the  latter  a  supersonic  wedge  flow.  This  report 
presents  results  using  that  same  combination  of  smoothing  techniques  applied 
to  flows  where  multiole  discontinuities  arise.  As  in  Reference  l.^lfhe  full 
inviscid  equations  of  motion  (Euler  equations),  cast  in  conservation  law 

form,  are  used  together  with  an  Adams -Bashforth  time  differencing  algorithm. 

'  n  -  f 

Two  classes  of^ multiple  discontinuity  inviscid  flows  are  solved:  (1)  a 
bursting  diaphragm  problem,  in  which  a  shock  wave  and  contact  surface  dis¬ 
continuity  are  simultaneously  present,  but  neither  have  yet  reached  a 
boundary,  and  (2)  the  flowfield  which  arises  when  two  normal  shock  waves  of 
unequal  strengths,  traveling  towards  each  other,  collide  and  give  rise  to 
two  shock  waves  of  new  and  different  strengths  along  with  a  contact  surface 
discontinuity. 

II.  Solution  Technique 

The  time -dependent,  one-dimensional,  Euler  equations  of  motion  in 
conservation  form  are  given  by: 


30  jl 
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and  p  is  the  pressure,  p  the  density,  u  the  velocity,  and  e  the  energy. 
Equation  (la)  is  solved  in  the  following  manner.  The  spatial  derivative  is 
obtained  pseudospectrally.  Starting  with  known  values  of  the  vector 
elements  of  E  at  a  specified  time,  a  Chebyshev  series  is  fitted  to  the 
values  of  each  E  element.  The  Chebyshev  series  is  represented  by 

N 

E(x,t)  =  Z  A  (t)T  (x )  ,*  (2) 

n  n  n 
n=0 

where  the  Chebyshev  functions,  Tn(x),  are  defined  by 


T^(x)  =  cos  [n  cos-1  (x)]  . 


(3) 


The  collocation  points  where  the  numerical  values  of  E  are  known  are  given 

by 


X. 
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0  <  j  <  N. 


(4) 


Therefore,  equation  (2)  becomes 

E(xj’°  “n^O  An(t)  COS  ^N1*  (5) 

To  fit  a  Chebyshev  series  to  the  E  data  requires  the  solution  of  (5)  for  the 
An's.  This  is  done  using  an  inverse  FFT. 

The  spatial  derivative  is  then  determined  from  the  ^'s  in  the 
following  way.  A  different  Chebyshev  series  is  used  to  represent  the 
spatial  derivative. 
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f—  *  E  A  (1)(t)  T  (x),  (6) 

3x  n  n 

n=0 

where 

(1)  2  N_1 

An  =  2  PAn*  (7) 

n  cn  p=n+l  P 

p+n*odd 

C  *  2,  C  *  1  for  0  <  n  <  N  , 

0  n 

0.0  . 

Performing  the  sum  in  (7)  for  the  A  ^'s  and  using  a  direct  FFT  in  (6) 

tl 

yields  the  spatial  derivative  values  at  each  of  the  collocation  points. 

The  solution  for  U  is  then  advanced  in  time  by  using  the  Adams 
Bashforth  algorithm. 

>  t-Hit  =  +  t  +  3  _  1  At(|Ej  +D  (8) 

1  j  2  k3x;  j  2  '•3x>j  j 

This  process  (equations  5  through  8)  is  then  cyclically  repeated.  The 
dissipation  term  Dj  is  evaluated  from  the  following  finite  difference 
representation  of  a  fourth  derivative. 

D,  ■  <"j+2+  °i-2-4  +  6  •• 

where  j  denotes  the  spatial  position  (collocation  point)  index.  The 
spectral  filter  (reference  4),  used  to  damp  out  the  high  frequency  solution 
components,  is  given  by 
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where  Kn  is  the  soectral  wave  number,  is  the  maximium  wave  number 

(corresponding  to  the  total  number  of  collocation  points)  and  K.  *  4  K 

o  6  max 
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The  time  integration  step  size  was  determined  from  the  pseudospectral 
form  of  the  Courant-Fredrich-Lewy  (CFL)  condition,  with  a  Courant  number  of 

0.5. 
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Here  u  is  the  velocity,  N  is  the  total  number  of  points,  c  is  the  speed  of 
sound,  and  CN  is  the  Courant  number. 


The  boundary  conditions  were  based  on  the  flow  which  arose  at  each 

boundary.  For  the  bursting  diaphragm  problem,  there  was  no  flow  (velocity  = 

zero)  at  the  boundaries.  (The  calculations  were  terminated  before  the  shock 

or  expansion  fan  reached  a  boundary.)  Therefore,  the  boundary  conditions 

were  to  keep  all  physical  variables  held  fixed  throughout  the  computation. 

For  the  colliding  shock  waves  problem,  variables  were  held  fixed  at 

supersonic  inflow  boundaries.  This  is  a  physical  boundary  condition  since 

the  integration  times  were  kept  well  below  those  where  either  the  contact 

surface  or  shock  wave  reached  their  respective  boundaries.  Outflow  boundary 

conditions,  namely  the  derivative  equal  to  zero,  were  not  considered  in  the 

present  work.  However,  at  subsonic  inflow  boundaries,  a  different  approach 

was  used.  At  a  subsonic  inflow  boundary  two  inflow  characteristics  exist, 

namely  ^  =  u  and  u  +  c.  The  single  outgoing  characteristic  propagates 
,dt 

along  *  u-c.  The  respective  characteristic  values  are  (p  -pc2),  (p  +  puc) 

and  (p  -  puc),  where  c  is  the  speed  of  sound.  Following  reference  5  two  flow 
variables  are  specified  at  the  subsonic  inflow  boundary,  and  the  third  is 
computed  from  the  calculated  value  of  the  outgoing  characteristic. 


III.  BURSTING  DIAPHRAGM  RESULTS 

Two  cases  were  considered: 

Density  Ratio 
2  to  1 
8  to  1 


Pressure  Ratio 

(1)  5  to  1 

(2)  10  to  1 
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They  will  be  discussed  sequentially . 

Results  for  case  (1)  are  shown  in  figures  1  through  4.  The  diaphragm 
separating  the  two  stagnation  zones  is  located  at  x  =  0.0.  One  hundred 
twenty-eight  Chebyshev  series  terms  were  used  to  represent  the  flow.  The 
time  integration  was  carried  out  to  t  =  0.378  (3500  time  steps).  The 
analytic  solution  is  reoresented  by  the  solid  line.  The  dashed  line 
represents  the  analytic  position  of  the  contact  surface  discontinuity.  The 
angled  solid  line  represents  the  expansion  fan.  At  all  times  the  shock  wave 
is  correctly  represented  as  a  discontinuity  traveling  with  the  proper 
velocity  (no  phase  error).  The  contact  discontinuity  is  spread  uniformly 
over  four  to  five  grid  points.  That  is,  the  points  all  fall  on  a  straight 
line  whose  slope  is  slightly  less  than  ninety  degrees  (the  exact  discon¬ 
tinuity  value).  The  same  behavior  is  present  in  the  energy  plots.  It  is 
not  as  clearly  visible  in  the  figures  because  the  difference  in  energies  on 
the  two  sides  of  the  contact  surface  front  are  very  small  in  comparison  to 
the  difference  in  densities.  The  pressure  and  velocity  correctly  show  no 
variation  across  the  contact  front.  The  expansion  fan  is  represented 
properly.  Some  minor  humps  are  present  both  on  the  high  and  low  state 
variable  sides,  and  there  is  a  slight  mismatch  in  slope.  This  mismatch  is  a 
function  of  the  physical  variable  being  plotted.  It  is  greatest  for  the 
velocity  while  least  (and  about  the  same)  for  the  remaining  three  flow 
variables. 

Results  for  case  2  are  shown  in  figures  5  through  8.  The  attributes 
discussed  above  are  also  present  here,  as  is  expected.  In  the  velocity 
plots  there  are  two  zones  of  oscillations  in  the  region  behind  the  shock 
wave.  It  turns  out  that  the  dividing  line  is  the  contact  surface 
location.  In  front  of  the  contact  surface  these  oscillations  are  very 
small,  while  behind  it  they  are  non-existent.  Additional  dissipation  would 
remove  this  imperfection  but  at  the  expense  of  a  minor  loss  of  shock 
resolution. 


IV.  COLLIDING  SHOCK  HAVES  PROBLEM 


To  Che  authors's  knowledge  this  is  the  first  time  this  problem  has  been 
treated  by  pseudospectral  means.  Two  cases  are  considered: 

Left  Hand  Side  Shock  Right  Hand  Side  Shock 

(1)  M  =  2.5  M  -  1.5 

(2)  4.0  1.5 

The  right  hand  side  shock  wave  is  propagating  to  the  left  while  the  left 
hand  side  shock  wave  propagates  to  the  right.  These  two  cases  were  chosen 
because  after  the  shock  waves  intersect,  the  shock  on  the  left  (initially  on 
the  right)  moves  in  different  directions  with  respect  to  ground  fixed 
coordinates  for  each  case.  It  moves  to  the  left  for  case  1  and  to  the  right 
for  case  2.  These  cases  result  in  supersonic  inflow  at  the  left  hand  side 
computational  boundary  and  subsonic  inflow  at  the  right  hand  side  computa¬ 
tional  boundary.  The  boundary  conditions  were  to  keep  all  flow  variables 
fixed  at  supersonic  inflow  points,  while  at  subsonic  inflow  points  the 
Dressure  and  velocity  were  specified  and  the  characteristic  value  of  the 
incoming  characteristic  was  used  to  calculate  the  density.  The  energy 
equation  is  used  to  calculate  the  energy.  (Even  though  characteristic 
boundary  conditions  were  used  for  the  colliding  shock  problem,  this  fourth- 
order  smoothing  was  not  able  to  control  Gibbs  oscillations  emanating  from 
the  subsonic  inflow  boundary.  A  second-order  scheme  was  required  in  the 
neighborhood  of  that  boundary  to  keep  those  oscillations  under  control. 
Without  it  the  computed  solution  rapidly  diverged  from  the  analytic  solution 
at  the  boundary.) 

The  initial  conditions  are  shown  in  Figures  9  through  12.  Results  for 
the  first  case  are  shown  in  Figures  13  through  24  .  The  analytic  solution 
is  shown  for  comparison  by  the  solid  lines.  For  this  case  the  shock  waves 
intersect  at  t  »  .057  .  These  conditions  have  been  chosen  so  that,  after 
the  intersection,  the  shock  on  the  left  is  propagating  to  the  left  with 


respect  to  ground-fixed  coordinates  in  case  1  and  to  the  right  in  case  2. 

The  arrows  which  appear  in  ail  post-intersection  figures  are  used  only  to 
indicate  the  direction  of  propagation  of  each  shock  with  respect  to  a  ground 
fixed  frame  of  reference  and  bear  no  relation  to  other  shock  character¬ 
istics.  As  previously  mentioned  the  two  post-intersection  shock  waves 
travel  in  opposite  directions  for  this  case. 

The  solution  at  post  collision  times  is  plotted  in  Figures  17  through 
24,  corresponding  to  the  2000th,  4000th  and  6000th  iterations  respectively. 
The  solid  lines  represent  the  analytic  solution.  Agreement  in  the  position 
of  the  discontinuities  between  the  analytic  and  computed  solutions  is  very 
good.  At  the  2000th  iteration  the  agreement  with  the  flow  variables  in  the 
region  between  the  two  shock  waves  is  only  fair  since  only  five  points  lie 
in  this  zone  at  this  time.  However,  as  the  number  of  iterations  increase 
and  this  zone  becomes  larger  (i.e.  more  points  lie  within  it),  the  agreement 
becomes  very  good.  By  the  4000th  iteration,  the  computed  values  of  flow 
variables  exactly  match  the  analytic  values. 

The  difference  in  resolution  of  the  two  shock  waves  is  due  to  two 
factors.  The  first  is  that  the  artificial  viscosity  constant  cannot  be 
tailored  to  both  shocks,  which  are  of  different  strengths.  The  value  used 
herein  was  selected  to  provide  maximum  resolution  of  the  stronger  shock 
while  producing  minimum  rounding  of  the  weaker  shock.  As  can  be  seen  in  the 
figures  there  is  some  minor  rounding  of  the  weaker  shock  wave.  The  second 
and  more  significant  reason  is  that  the  physical  space  point  resolution  is 
not  constant.  It  is  greatest  at  x  =  ±  1  and  least  at  x  =  0.0.  The  weaker 
shock  in  this  case  always  lies  in  the  neighborhood  of  x  =  0.0,  the  minimum 
point  resolution  zone.  This  is  why  the  computed  shock  front  of  the  weaker 
shock  is  not  as  steep  as  for  the  stronger  shock  wave. 

The  oseudospectral  resolution  characteristics  of  the  contact  discon¬ 
tinuity  are  not  as  good  as  those  of  the  shock  waves.  Part  of  this  is  due  to 
the  above  mentioned  unequal  point  spacing.  However,  most  is  due  to  the 
nature  of  the  solution  scheme  together  with  the  significant  difference  in 
the  physical  nature  of  the  two  types  of  discontinuities.  Shock  waves  yield 
discontinuous  values  of  all  flow  variables  across  their  front  while  contact 


surface  discontinuities  yield  jumps  in  density  and  energy  alone. 

Results  for  the  second  case  are  shown  in  Figures  25  through  40  .  There 
are  several  important  differences  between  this  case  and  the  first.  The 
initial  shock  strengths  are  such  that,  after  intersection,  both  shock  waves 
move  to  the  right  with  respect  to  ground  fixed  coordinates.  Further,  the 
difference  in  post  intersection  shock  strengths  is  greater  in  this  case. 
Perhaps  the  major  difference,  however,  is  in  the  initial  positions  of  the 
shocks.  In  the  first  case  both  were  far  away  from  the  computational 
boundaries.  In  this  case  it  was  decided  to  purposely  place  one  shock  much 
closer  to  a  boundary  to  see  what,  if  any,  effect  this  would  have  on  the 
properties  of  the  solution.  The  left  hand  shock  was  placed  only  twenty 
points  from  the  computational  boundary.  No  adverse  results  were  observed. 

As  in  the  first  case,  the  comparison  with  the  analytic  solution  is  very 
good.  All  shocks  are  resolved  as  sharp  discontinuities  and  at  the  correct 
position.  The  contact  surface  resolution  was  as  before  not  as  good,  being 
spread  over  four  to  five  collocation  points,  though  centered  about  the 
analytic  location. 

V.  CONCLUSIONS 

The  central  conclusion  is  that  the  pseudospectral  solution  technique 
together  with  a  fourth-order  artificial  viscosity  smoothing  scheme  works 
just  as  well  for  cases  where  multiple  discontinuities  are  present  as  it  does 
when  only  a  single  discontinuity  is  present.  The  resolution  of  contact 
surface  discontinuities  is  not  as  good  as  that  for  shock  discontinuities, 
the  former  being  spread  over  four  to  five  grid  or  collocation  points  as 
compared  to  two  to  three  grid  points  (one  to  two  grid  intervals  for  the 
latter ) . 
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Figure  1  Bursting  Diaphragm  Flow  Pressure  Field,  Case 


Figure  2  Bursting  Diaphragm  Flow  Density  Field,  Case  i. 


Figure  3  Bursting  Diaphragm  Flow  Velocity  Field,  Case 
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Figure  4  Bursting  Diaphragm  Flow  Energy  Field,  Case 


13 


Figure  5  Bursting  Diaphragm  Flow  Pressure  Field,  Case 
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Figure  7  Bursting  Diaphragm  Flow  Velocity  Field,  Case 
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Figure  8  Bursting  Diaphragm  Flow  Energy  Field,  Case 
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Colliding  Shock  Waves,  Initial  Density  Field,  Case 
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Figure  12  Colliding  Shock.  Waves,  Initial  Energy  Field,  Case 


POST- INTERSECT  I ON 


Figure  14  Colliding  Shock  Waves,  Post-Collision  Density  Field,  Case 
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Figure  15  Colliding  Shock  Waves,  Post-Collision  Velocity  Field,  Case 
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Figure  18  Colliding  Shock  Waves,  Post-Collision  Density  Field,  Case 


Figure  19  Colliding  Shock  Waves,  Post -Collision  Velocity  Field,  Case 
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Colliding  Shock  Waves,  Post-Collision  Energy  Field,  Case 
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Figure  21  Colliding  Shock  Waves,  Post-Collision  Pressure  Field,  Case 
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Figure  24  Colliding  Shock  Waves,  Post-Collision  Energy  Field,  Case 
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Colliding  Shock  Waves;  Initial  Pressure  Field,  Case 
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Figure  26  Colliding  Shock  Waves;  Initial  Density  Field,  Case 


gure  28  Colliding  Shock  Waves;  Initial  Energy  Field,  Case 


figure  29  Colliding  Shock  Waves;  Post-Collision  Pressure  Field,  Case 
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Figure  31  Colliding  Shock  Waves;  Post-Collision  Velocity  Field,  Case 
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figure  32  Colliding  Shock  Waves;  Post-Collision  Energy  Field,  Case 
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<’igure  33  Colliding  Shock  Waves;  Post-Collision  Pressure  field,  Case 
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Figure  36  Colliding  Shock  Waves;  Post-Collision  Energy  Field,  Case 
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■’igtire  37  Colliding  Shock  Waves;  Post -Collision  Pressure  Field,  Case 
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Figure  38  Colliding  Shock  Waves;  Post-Collision  Density  Field,  Case 
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Figure  40  Colliding  Shock  Waves;  Post-Collision  Energy  Field,  Case 
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